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Abstract 

Background: Degradation of cellulose to glucose requires the cooperative action of three classes of enzymes, 
collectively known as cellulases. Endoglucanases randomly bind to cellulose surfaces and generate new chain ends 
by hydrolyzing (3-1,4-D-glycosidic bonds. Exoglucanases bind to free chain ends and hydrolyze glycosidic bonds in 
a processive manner releasing cellobiose units. Then, (3-glucosidases hydrolyze soluble cellobiose to glucose. 
Optimal synergistic action of these enzymes is essential for efficient digestion of cellulose. Experiments show that as 
hydrolysis proceeds and the cellulose substrate becomes more heterogeneous, the overall degradation slows down. 
As catalysis occurs on the surface of crystalline cellulose, several factors affect the overall hydrolysis. Therefore, 
spatial models of cellulose degradation must capture effects such as enzyme crowding and surface heterogeneity, 
which have been shown to lead to a reduction in hydrolysis rates. 

Results: We present a coarse-grained stochastic model for capturing the key events associated with the enzymatic 
degradation of cellulose at the mesoscopic level. This functional model accounts for the mobility and action of a 
single cellulase enzyme as well as the synergy of multiple endo- and exo-cellulases on a cellulose surface. The 
quantitative description of cellulose degradation is calculated on a spatial model by including free and bound 
states of both endo- and exo-cellulases with explicit reactive surface terms (e.g., hydrogen bond breaking, covalent 
bond cleavages) and corresponding reaction rates. The dynamical evolution of the system is simulated by including 
physical interactions between cellulases and cellulose. 

Conclusions: Our coarse-grained model reproduces the qualitative behavior of endoglucanases and exoglucanases 
by accounting for the spatial heterogeneity of the cellulose surface as well as other spatial factors such as enzyme 
crowding. Importantly, it captures the endo-exo synergism of cellulase enzyme cocktails. This model constitutes a 
critical step towards testing hypotheses and understanding approaches for maximizing synergy and substrate 
properties with a goal of cost effective enzymatic hydrolysis. 

Keywords: Cellulose degradation, Synergy, Exo-cellulase, Endo-cellulase, Agent-based model, Spatial heterogeneity 



Background 

Biofuel production from lignocellulosic materials is con- 
sidered to be a promising option to substantially reduce 
the dependence on petroleum [1-3]. The conversion 
of lignocellulosic biomass (agronomic residues, paper 
wastes, energy crops) into ethanol consists of the extrac- 
tion and pretreatment of cellulose from the biomass, 
hydrolysis (the enzymatic breakdown of crystalline cel- 
lulose fibers into monomer glucose) and finally the 
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fermentation of glucose to ethanol. Current approaches 
mainly differ from one another in the method of pre- 
treatment. Cost-competitive production of ethanol is 
currently prevented by the low efficiency of converting 
cellulose into glucose [4]. Greater efficiency may be 
achievable through improvements in hydrolysis. 

Enzymatic hydrolysis of cellulose is a complex reac- 
tion. In the classical model, the heterogeneous catalytic 
cleavage of the glycosidic bond takes place on a crystal- 
line cellulose surface and requires the cooperative action 
of three classes of aqueous enzymes, collectively known 
as cellulases. These are (i) endoglucanases, (ii) exogluca- 
nases or cellobiohydrolases and (iii) p-glucosidases. 
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Recently, it has been proposed that oxidative enzymes 
(monoxygenases) may also play a role in cleaving glyco- 
side bonds, although this new mechanism may be 
restricted to certain types of microbes [5], It is widely 
accepted that endoglucanases cleave [3-1,4-D-glycosidic 
bonds at random sites within both amorphous and crys- 
talline polysaccharide chains, creating new chain ends 
on the cellulose surface [6-11]. Exoglucanases prefer to 
hydrolyze crystalline cellulose chains by acting on the 
free chain ends and releasing cellobiose units in a proces- 
sive manner. Soluble cellobiose units are then converted 
into glucose by p-glucosidases. Consequently, these 
enzymes display strong synergy [12-18]. 

The classical chemical kinetics assumption of uni- 
formly mixed systems does not hold in the case 
of enzymatic hydrolysis of cellulose fiber, as it is hetero- 
geneous in nature. Such reactions are rather character- 
ized by time-dependent rate constants and non-uniform 
concentration variation of reacting species. Although 
kinetic models [11,19-22] have been used to explain vari- 
ous features of the enzymatic hydrolysis of cellulose, they 
fail to account for spatial details of the cellulose substrate 
as well as the specificity of binding sites. Recently, 
Zhou and colleagues [23-25] proposed a "morphology- 
plus-kinetics rate equation approach" that explicitly 
captures the hydrolytic evolution of cellulose substrate. 
In addition, a kinetic model was developed based on 
population-balance equations in which a distribution of 
chains with different chain-lengths was explored [26,27]. 
Yet, these models give little insight regarding the action 
of cellulases at the molecular level. 

It is imperative to develop spatial models of cellulose 
degradation because spatial effects such as enzyme 
crowding on the cellulose surface have been shown to 
lead to a reduction in hydrolysis rates. In order to ac- 
count for the spatial heterogeneity of the system during 
cellulose hydrolysis, a cellular automata model [28] was 
developed to study the effect of different parameters such 
as enzyme binding and hydrolysis on the overall kinetics 
of cellulose by the cellulases. Alternatively, all-atom mo- 
lecular dynamics (MD) simulations can provide details 
of molecular level events at high precision. Recent 
MD simulation studies [29-32] have proven effective 
for understanding enzyme-substrate binding, proces- 
sivity and activity. However, because of length and 
time scale limitations, it is not currently possible to 
simulate the entire crystalline cellulose degradation 
process using all-atom MD simulations. 

We have developed a coarse-grained stochastic model 
that captures the interaction of endo- and exo-cellulases 
with crystalline cellulose at a mesoscopic level. This 
model was specifically designed to improve our under- 
standing of the molecular-level details of the enzymatic 
hydrolysis of crystalline cellulose. This paper introduces 



the basic framework and demonstrates how this model 
can be an effective and easily modifiable testing platform 
for new hypotheses based on experimental data on vari- 
ous cellulase components and substrate characteristics. 
By capturing the reactive nature of the cellulose sub- 
strate and the activities of non-complexed cellulases at 
the molecular level, this method forms a bridge between 
all-atom MD studies and deterministic reaction-rate 
approaches. To the best of our knowledge, it is the first 
model that is able to relate the synergetic action of mul- 
tiple enzymes to molecular level details such as the 
hydrogen bond network of a cellulose substrate. 

Results and discussion 

Model development 

The overall efficiency of the heterogeneous catalysis that 
occurs in the enzymatic hydrolysis of crystalline cellu- 
lose depends on factors such as adsorption, desorption, 
diffusion rates on the insoluble cellulose substrate, and 
processivity. In our model, catalysis is broken down into 
distinct parts related to different kinetic events (chemical 
reactions) performed by individual particles (enzymes). 
Specifically, we include the following reactive events: ad- 
sorption of cellulases on the solid cellulose substrate, 
inter-chain hydrogen bond breaking, hydrolysis of glyco- 
sidic bonds, and desorption of cellulases from cellulose. 
These reactions constitute the main elements of this 
model, and their realization is achieved by following and 
updating the state (based on certain predefined rules dis- 
cussed below) of each individual particle in the system 
as it evolves in time. The actions of cellulases are mod- 
eled based on the most abundant endoglucanase (EG I) 
and the two cellobiohydrolases (CBH I and CBH II) 
secreted by the filamentous fungus Trichoderma reesei 
[33], as these three enzymes have been widely studied 
[10,11,34-36] and have been the target for improvement/ 
design for efficient biodegradation [37,38]. 

Model of cellulose substrate 

The cellulose surface layer (Figure 1) is modeled as a two- 
dimensional grid, consisting of multiple glucan chains, 
each of them having the same number of monomers, i.e., 
the degree of polymerization of glucan chains is the same. 
Glucose molecules are linked to each other through intra- 
chain covalent glycosidic bonds, while links between inter- 
chain glucose molecules correspond to multiple hydrogen 
bonds. In this representation, all chains are oriented such 
that the left side corresponds to the nonreducing end and 
the right side to the reducing end of the glucan chains. 
This chosen orientation represents the most commonly 
found form of crystalline cellulose in nature, cellulose 1(3 
[39]. The state of each glucose unit is defined by a set of 
seven binary parameters, P lt P2, ■ ■ - J?7, listed in Table 1. 
Each row corresponds to one of the seven parameters 
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Glycosidic bond 



Hydrogen bond 

Figure 1 Model for the cellulose surface composed of glucose. 

At each time step, the stochastic simulation stores the intra- and 
inter-chain neighbors of each glucose unit. Additionally, each unit is 
characterized by a set of P rn parameters, 1 <m<7, listed in Table 1 . 



characterizing one glucose unit. The values of the binary 
parameters (0 and 1) represent various conditions of one 
glucose unit. 

The first parameter, Pj informs whether the glucose 
unit belongs to the cellulose substrate ('nonsoluble') or 
is in the aqueous phase ('soluble'). Initially, all glucose 
molecules are part of the cellulose surface (Pj=0), and 
as the simulation progresses, they become soluble 
{Pi = 1) in the form of simple sugars: glucose, cellobiose 
or cellotriose. The length of the soluble oligomers can 
be easily modified. 

The second parameter, P 2 specifies whether the glu- 
cose unit constitutes the nonreducing end (NE) or the 
reducing end (RE) of a glucan chain. Also, when a glyco- 
sidic bond in the middle of a chain is hydrolyzed, two 
new ends are created, one new NE and one new RE. 

Parameter P 3 informs whether an endo-cellulase cov- 
ers the glucose unit. Similarly, the values of parameters 
P 4 and P s indicate whether an exo-cellulase hydrolyzing 
from the reducing end (exo-R) or an exo-cellulase 
hydrolyzing from the nonreducing end (exo-N) covers 
the glucose unit. At any time, only one cellulase is 



Table 1 State variables of a glucose unit 3 


P/Value 


0 


1 


Pi 


nonsoluble 


soluble 


?2 


N.E. 


R.E. 


P 3 


uncovered 


covered by endo 


Pa 


uncovered 


covered by exo-R 


P, 


uncovered 


covered by exo-N 


?e 


not locked 


locked by exo-R 


P? 


not locked 


locked by exo-N 



allowed to cover a specific glucose, during which glyco- 
sidic and hydrogen bonds may be cleaved or broken. A 
glucose unit, which is not covered by any cellulases, may 
become locked by a processive cellulase, exo-R or exo-N. 
This is specified by parameters Pg and P7, respectively. A 
locked glucose unit only facilitates the binding of the 
locking exo-R (exo-N); it does not constitute an available 
binding site for any other cellulases, nor can its glyco- 
sidic or hydrogen bonds be cleaved or broken, until the 
locking cellulase binds directly to it. 

Model of cellulase with endo-activity 

Cellulases with endo-activity (referred to as endo- 
cellulases or simply endo) are modeled through a set of 
interactions between the cellulose surface and among 
themselves. A detailed description of the actions of 
adsorbed endo-cellulases is presented in Figure 2, while 
Table 2 lists the parameters that determine their overall 
activity. Additionally, a state parameter is used to specify 
whether the cellulase is adsorbed to the substrate or is 
in solution. In the future, we plan to incorporate another 
state parameter to describe a decrystallization step that 
"prepares" the substrate for productive binding. 

Endo-cellulases, once adsorbed to the cellulose surface, 
can break inter-chain hydrogen bonds, hydrolyze glyco- 
sidic bonds and desorb from the substrate into solution. 
Each of these chemical reactions is essentially a Poisson 
process that takes place at a specific, constant rate defined 
by the propensity function of that reaction [40] . 

Each endo-cellulase adsorbs to a randomly chosen 
available site. A site is a set of glucose units, as pre- 
sented in Figure 3a, consisting of nine consecutive glu- 
cose units in three neighboring chains, for a total of 27 
adjacent glucose units [35]. However, we assume that 
just a length of 4 glucose units in the middle chain is 
enough for it to form a productive complex. This choice 
was motivated by the empirical studies of Claeyssens 
et al. [41] and Biely et al. [42] who argued that the sub- 
strate binding site of EG I is an extended one, consisting 
of four sugar binding subsites with a catalytic group 
located in the middle. A site is available if all monomers 
in the middle glucan chain belong to the cellulose sub- 
strate and none of the twelve monomers is covered by 
another enzyme nor are they locked. Desorption of the 
cellulase might take place at any time. If the glycosidic 
bond in the catalytic region is already hydrolyzed, the 
cellulase desorbs into solution within an exponentially 
distributed time interval with rate parameter k 0 ^ ast (see 
Figure 2). The catalytic region of an endo-cellulase is 
considered to be the glycosidic bond between the second 
and third glucose units in the middle chain covered by 
the enzyme, shown in Figure 3c. The hydrolysis of the 
glycosidic bond can only take place after all inter-chain 
hydrogen bonds between the covered glucose units are 
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Figure 2 Flowchart presenting the actions of an endo-cellulase. The actions are: adsorption onto the cellulose surface (ADSORB), breaking 
hydrogen bonds (BREAK H. Bs.), hydrolyzing the glycosidic bond (CLEAVE Glyc. BJ, and desorption from the cellulose (DESORB). Green rectangles 
denote chemical reactions (events) and orange ellipses denote branching points. Abbreviations: 'Glyc. B.' covalent glycosidic bond; H Bs.' 
hydrogen bonds; 'nrHB' number of hydrogen bonds present between the monomers covered by the cellulase. 



broken (Figure 3b). The time it takes to break the hydro- 
gen bonds is proportional to the number of bonds that 
need to be broken, denoted by 'nrHB' in Figure 2. If 
there is at least one hydrogen bond that needs to be 
broken, the cellulase breaks it or desorbs into solution. 
Similarly, if all hydrogen bonds have already been 
broken, as shown in Figure 3b, the cellulase either 
hydrolyzes the glycosidic bond at the catalytic region or 
desorbs into solution. These decisions are implemented 
using Gillespie's algorithm [40]. Hydrolysis of a glyco- 
sidic bond is always followed by desorption of the cellu- 
lase within an exponentially distributed time with rate 
parameter k 0 g. 

Models for celluloses with exo-activity 

Cellulases with exo-activity can be of two types: cellu- 
lases hydrolyzing the glucan chain from its reducing end 
(referred to as exo-R cellulase or simply exo-R), and 
cellulases hydrolyzing the glucan chain from its non- 
reducing end (referred to as exo-N cellulase or simply 
exo-N). The interactions between adsorbed exo- 

Table 2 Rate constants characterizing endo-cellulases 

Nomenclature 



k on adsorption rate constant 

k olf desorption rate constant 

k of ff ast a higher desorption rate constant than k off 

khbbreak rate constant for breaking a single hydrogen bond 

k g iy rate constant for hydrolyzing a glycosidic bond 



cellulases and the cellulose surface are shown in 
Figure 4, while Table 3 lists the relevant parameters 
that determine the overall activity of exo-cellulases. 
Additionally, a state parameter is used to specify 
whether the cellulase is adsorbed to the substrate or is 
in solution. Again, we intend to consider an additional 
state parameter to describe the decrystallization step in 
future models. 

Exo-cellulases, once adsorbed to the cellulose surface, 
can break hydrogen bonds, hydrolyze glycosidic bonds, 
slide along a chain, and desorb from the cellulose into 
solution. For simplicity, in the following description we 
restrict ourselves to actions carried out by an exo-R cel- 
lulase. The ones carried out by an exo-N cellulase are 
essentially the same. 

The adsorption site of an exo-R cellulase consists of 
nine consecutive glucose units in three neighboring 
chains, for a total of 27 adjacent glucose units (shown in 
Figure 5a), subject to the condition that the middle glu- 
can chain has a reducing end. This choice has been 
motivated by the three-dimensional structure of the 
catalytic domain of CBH I from T. reesei [34,43]. This 
catalytic site resides within a relatively long (~50 A) cel- 
lulose binding tunnel holding ten glucose molecules, out 
of which three — near the outlet — form the product bind- 
ing sites. As cellobiose is the main product released by 
CBH I [16,44], the exo-cellulase in our model has only 
two product binding sites, which, for the substrate, 
translates into a total of nine sugar binding sites. The 
adsorption site for an exo-cellulase is considered to be 



Asztalos et al. Biotechnology for Biofuels 2012, 5:55 
http://www.biotechnologyforbiofuels.eom/content/5/1/55 



Page 5 of 1 5 



Adsorption 



X Broken hydrogen bond 
X Cleaved glycosidic bond 




Complexation 




Cleaving glycosidic bond 




Figure 3 Interactions between endo-cellulase and cellulose crystals. Schematic representation of (a) an endo-cellulase adsorbed onto the 
cellulose surface, (b) hydrogen bonds breaking between the monomers covered by the enzyme and (c) hydrolysis of the glycosidic bond. 



available if all monomers in the middle glucan chain be- 
long to the cellulose substrate and none of the 27 mono- 
mers is covered by another enzyme nor are they locked. 

During adsorption, the inter-chain hydrogen bonds of 
the glucose units beneath the cellulase are instantly 
broken (Figure 5a). Similar to an endo-cellulase, desorp- 
tion of an exo-cellulase from the cellulose sheet can take 
place at any time. If any of the glycosidic bonds between 




[STAY | [PROCESS] [ DRSORB | [ DESORB ] [ DESORB | 
I 



Figure 4 Flowchart presenting the actions of an exo-cellulase 
after adsorption to the cellulose surface. Chemical reactions are 
denoted by green rectangles; orange ellipses denote branching 
points. Abbreviations are the same as used in Figure 2. Additionally, 
r is a uniformly, distributed random variable between 0 and /. 



the glucose units covered by the cellulase in the middle 
glucan chain is already hydrolyzed, the cellulase desorbs 
from the surface within an exponentially distributed 
time interval with rate parameter k 0 jff ast . This assump- 
tion was built into the model in order to account for the 
continuity of the chain entering the tunnel of CBH I. 
The frequency of cellulase dissociations from the sub- 
strate is set by a probability a. Results of high speed 
AFM measurements [45] indicate that once the cellulase 
is adsorbed to a free end, it continues to process the 
chain until it reaches the end of the chain. In this light, 



Table 3 Rate constants and relevant parameters for 
characterizing and exo-cellulase 

Nomenclature 

k 0 ff 

kafffast 

a 



thbbreak 



adsorption rate constant 

desorption rate constant 

a higher desorption rate constant than k off 

probability for the exo-cellulase to desorb 
from cellulose 

time for an exo-cellulase to break a single 
inter-chain hydrogen bond 

time for an exo-cellulase to hydrolyze a 
glycosidic bond and slide one cellobiose 
unit along the glucan chain 

time for an exo-cellulase to remain at a 
certain location 
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X Broken hydrogen bond 
X Cleaved glycosidic bond 



Adsorption and Complexation 




Processive degradation of one glucan chain 




Figure 5 Interactions between exo-cellulase and cellulose crystal. Schematic representation of (a) an exo-R cellulase adsorbed to the 
cellulose surface followed immediately by the breaking of hydrogen bonds between the monomers covered by the cellulase and (b) the 
processivity of the glucan chain by an exo-R cellulase; it comprises the hydrolysis of the glycosidic bond and cellulase directed movement along 
the chain. 



a is usually set to a relatively low value. The 'Spatial 
Conditions' branching point in Figure 4 checks whether 
the cellulase can further hydrolyze the glucan chain. 
When another cellulase obstructs its diffusion, it may 
stay in place for a time of length t stay . It may also desorb 
within an exponentially distributed time with parameter 
k 0 ff, if glucose units are missing from the middle chain 
or if it has reached the nonreducing end of the chain 
(and similarly for the exo-N). 

The processive movement of exo-cellulases is modeled 
as a two-step process: (i) the cellulase hydrolyzes the 
glycosidic bond positioned at the active site, followed by 
(ii) its sliding along the processed chain by one cello- 
biose unit, instantly breaking inter-chain hydrogen 
bonds [34] (Figure 5b). The active site of an exo-R is 
considered to be the glycosidic bond between the second 
and third glucose unit from the reducing end of the 
middle, processed glucan chain. If the two glucose units 
in front of the cellulase belong to the substrate, and 
none of the six glucose units in front of the cellulase 
(two consecutive ones in three chains) are occupied by 
any other cellulases, these six monomers become locked, 
implying that they could be covered by the exo-cellulase. 
Locked glucose units are not considered available bind- 
ing sites for endo- or exo-cellulases. The processing 
time of exo-cellulases is calculated as t process = t move + 
thbbreak'nrHb, where t move specifies the time during 
which the cellulase hydrolyzes the glycosidic bond at the 
active site and slides along the middle glucan chain by 
one cellobiose unit. Here nrHb is the number of present 
hydrogen bonds between the six locked monomers. The 
same strategy is employed for exo-N except that the 
processivity is along the opposite direction towards the 



left-hand side. Furthermore, an exo-cellulase is never 
allowed to productively bind to chains in the middle of 
the cellulose surface. The exo-R cellulases are only 
allowed to productively bind to a free reducing end, 



INITIALIZATION 
[E] 0 , Number of glucose units 
Kinetic Parameters, Time=0.0 



[ DECIDE THE FIRST EVENT ] 

[ UPDATE THE TIME*] 
1 

[ PERFORM THE EVENT V- 



DECIDE THE NEW EVENT ^ 
FOR THAT ENZYME 
Based on the flowcharts for endo & exo J 



CHECK FOR 
ADSORPTION 



ENTER THE EVENT IN THE 
MASTER QUEUE OF EVENTS J 



CHOOSE THE NEXT EVENT 
UPDATE THE TIME 



1 



J 



DEGRADATION TRESHOLD - 



jyes 

STOP 

Figure 6 Simulation timeline. Here an 'EVENT may refer to any 
chemical reaction that involves a cellulase (adsorption, desorption) 
or is catalyzed by a cellulase (e.g., inter-chain hydrogen bond 
breaking, hydrolysis of glycosidic bond). 
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while the exo-N cellulases are only allowed to product- 
ively bind to a free non-reducing end. 

Algorithm for time evolution 

The outline of the overall simulation is sketched in 
Figure 6. Time is measured in seconds and is advanced 
in a continuous and asynchronous manner. The simula- 
tion starts with the adsorption of an individual cellulase 
(endo or exo) and it proceeds following the flowcharts 
presented in Figure 2 (for endo-cellulase) or Figure 4 
(for exo-cellulase). Each individually adsorbed cellulase 
is followed separately over the course of the simulation. 
While in solution, they all compete with each other for 
adsorption, and once adsorbed, the selection of the 
chemical reactions they are involved in follows a well- 
defined, rule-based schematic (see below) that finally 
leads to the hydrolysis of the entire cellulose substrate. 
However, Gillespie's algorithm [40,46,47] plays a crucial 
part in the simulation, as many chemical reactions, once 
selected by the rule-based scheme, are modeled as Pois- 
son processes and therefore are implemented using this 
algorithm. Gillespie's algorithm is a Monte Carlo tech- 
nique that allows one to sample the ensemble of trajec- 
tories for a set of biochemical reactions. It models 
chemical reactions as a stochastic process and remains 
valid for low copy numbers of reactants. Here we imple- 
ment the direct version of this method [46] . 

Chemical reactions that involve endo-cellulases (ad- 
sorption, desorption) or are catalyzed by endo-cellulases 
(inter-chain hydrogen bond breaking, hydrolysis of 
glycosidic bond) are modeled as Poisson processes. This 
choice has been motivated by the knowledge of the spe- 
cific reaction rate constants. 

A different modeling approach is taken in the case of 
exo-cellulases. While adsorption of exo-cellulases to the 
cellulose substrate is modeled as a Poisson process, each 
adsorbing event is combined with an instantaneous 
structural transformation of the surface as hydrogen 
bonds below the cellulase are assumed to break at the 
very moment of adsorption. The processivity of exo- 
cellulases is defined through specific rules: when certain 
spatial conditions are met, the cellulase slides along the 
glucan chain by one cellobiose unit leaving behind a sol- 
uble cellobiose and instantly changing the cellulose sur- 
face (breaking intact hydrogen bonds) beneath itself. We 
assume that each exo-cellulase has the same, constant 
processing velocity so the only difference in processing 
times of one exo-cellulase from another originates in the 
number of intact hydrogen bonds covered by respective 
cellulases. Finally, desorption of exo-cellulases is also 
modeled as a Poisson process with a constant rate 
parameter. 

Each chemical reaction is modeled as a discrete event 
occurring instantaneously while the state of the system 



remains unchanged between two consecutive events. 
The events associated with each of the adsorbed cellu- 
lases are stored in a priority queue, here referred to as 
the master queue of events (Figure 6). They are sorted by 
the simulated time at which they should occur. The 
simulation runs until the surface degrades to a specific 
degradation threshold, which is usually set to 100%, or 
the point at which the substrate is not able to adsorb 
more enzyme particles. 

Simulation parameters 

The values of the kinetic parameters we employed are 
listed in Table 4. The values of k on and k a ff were deter- 
mined based on the adsorption equilibrium constants, 
which for both types of cellulases were equal to JO 3 M 1 . 
The fast desorption rate constant {k 0 jff asl ) was always one 
order of magnitude larger than k 0 jf. The values of the hy- 
drolysis rate constants for both endo- and exo-cellulases 
{kgiy) were the ones estimated by Zhang and Lynd [22] 
and accordingly, the rate constant for endo-cellulases 
was set at fivefold that used for exo-cellulases [11,48]. 
From there we obtain t move = 12 s; t stay was set to the 
same value. The probability a was set to 0.1 and the time 
to break a hydrogen bond to 10 12 s [49]. In the results 
presented in the next section, we did not include hydro- 
gen bond reformation, but this is a reaction that can eas- 
ily be included in this model. 

Table 4 also includes the molecular weights of EG I and 
CBH I used in the simulation. The molecular weights of 
glucose (ja = 180.15588 g/mol) and anhydroglucose 
(ja, = 162.1406 g/mol) molecules are essential in the calcula- 
tion of the soluble sugar concentration. The cellulose sub- 
strate is composed of five glucan chains, each of them 
having 4000-5000 glucose monomers. The area of a cel- 
lobiose unit [11] was set to A G2 = 5.512 x Iff 19 m 2 . The 
initial (molar) enzyme concentration was calculated from 
the number of cellulases present in the system and the 
volume of the system, V=Sd, where 5 is the cellulose sur- 
face area in m 2 and d is of the order of um. In most of 
the runs the initial enzyme concentration was set to be 
[E] 0 = 2 fiM. In the following section, unless stated other- 
wise, we used the parameter values listed in Table 4 and 
the above enumerated initial conditions. 



Table 4 Input parameters used in simulations unless 
otherwise stated 



Name 


Notation 


Endo 


Exo-R/Exo-N 


Molecular weight b 




52500 g/mol 


63500 g/mol 


Adsorption rate constant 


kon 


100 (sMf 1 


10 4 {sMr 1 


Desorption rate constant 


koff 


0.1 s" 1 


10 s"' 


Higher desorption rate constant 


kofffast 


1 s- 1 


1 00 s" 1 


Glycosidic bond hydrolysis c 


kgly 


0.35 s" 1 


0.0846 s" 1 
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The output consists of the time evolution of the con- 
centrations of glucose, cellobiose and cellotriose present 
in the aqueous phase, the adsorption density, and the 
number of available binding sites per gram of cellulose. 
Results were averaged over 10 replicas of the system. 
The conversion is quantified by counting the number of 
all the soluble glucose molecules including those in cel- 
lobiose and cellotriose and dividing that by the initial 
number of glucose molecules. 



Hydrolysis by endo-cellulases 

First we simulate and analyze the hydrolysis of a crystal- 
line cellulose layer solely by endo-cellulases. Figure 7a 
shows the timeline of percent cellulose degradation. 
After an initial slow hydrolysis phase, the simulation 
results agree qualitatively with published experimental 
results [50]. Using parameter values listed in Table 4, 
our model reproduces well the observed experimental 




Time(h) 

Figure 7 (a) Hydrolysis by endo-cellulases. Time course of 
hydrolysis by endo-cellulases.The inset shows the variation in the 
conversion time needed to degrade 5%, 50% and 80% of the 
substrate as a function of the adsorption rate constant k on . (b) Tota 
sugar production and individual monomer and oligomer 
component distribution over time. 



hydrolysis time scales. As expected, the time it takes to 
convert a given percent of the substrate decreases as k on 
increases (Figure 7a inset). Similarly, the relative produc- 
tion of soluble oligosaccharides shown in Figure 7b 
agrees well with the experiments [11]. Cellobiose is the 
major type of soluble sugar, both in experiments and in 
our model. In our model, the final molar ratio between 
glucose and cellobiose is close to 1:10, while the final 
molar ratio between cellotriose and cellobiose is close 
to 7:10. In experiments, the glucose concentration is 
observed to be higher than the cellotriose concentration, 
while our model shows the opposite. Our model assumes 
that cellulose oligomers of length less than 4 enter solu- 
tion and none of the modeled cellulases can digest them 
after they enter solution phase. This simplified assump- 
tion leads to disagreement with experiments. 

The effect of varying the initial [E] 0 cellulase concen- 
tration upon conversion times is plotted in Figure 8a. As 
the enzyme concentration increases, the gap between 




5 10 15 20 25 

Initial enzyme concentration (|J.M) 



30 




Time(h) 

Figure 8 (a) Effect of initial concentration of endo-cellulase on 
hydrolysis of cellulose. Change in conversion time as the initial 
enzyme concentration is varied. (b) Time course of enzyme 
adsorption. The inset shows the percentage of adsorbed enzymes as 
a function of time. The cellulose layer is composed of five glucan 
chains, each of them composed of 5000 glucose units. 
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the time to convert 50% and 75% of the cellulose sub- 
strate decreases. At high enzyme loading (> 22 piM), the 
conversion time reaches a constant value as the sub- 
strate is saturated by adsorbed enzymes: already bound 
enzymes mutually obstruct the adsorption of additional 
enzymes onto the surface. As the substrate is reduced 
over time, the number of available binding sites also 
decreases, thus fewer and fewer enzymes are able to 
bind to the surface. This trend is captured in Figure 8b. 
Cellulases quickly adsorb onto the surface at early times 
in the hydrolysis process (not shown because of the 
small time scale), after which their number follows a 
Poisson decay, and desorption is modeled as a Poisson 
process. The inset from Figure 8b shows the decrease in 
bound enzyme percentage as the initial enzyme concen- 
tration increases. Although this decrease is relatively 
small, it underlines the finite size effect of the substrate, 
which is consistent with Figure 8a. 

Although the simulation results show good agreement 
with the experiments, the lag phase observed in the hy- 
drolysis curve (Figure 7a) is unexpected, as it was not 
observed in any bulk measurements involving enzymatic 
hydrolysis of cellulose. It was, however, observed during 
acid hydrolysis of bacterial cellulose [51] and during an- 
aerobic bacterial digestion of cellulose [52]. The possible 
reasons behind the occurrence of the lag phase are vari- 
ous, (i) At the beginning, the random glycosidic bond 
cleavages are too far from each other to release glucose, 
cellobiose or cellotriose, as this requires a finite amount 
of time for cellulases to revisit the neighborhood of a 
cleaved bond, (ii) The model used here does not take 
into account enzyme diffusion [53], which could acceler- 
ate the ability of the enzyme to locate the neighborhood 
of a cleaved glycosidic bond. (Hi) In contrast to a realistic 
cellulose crystal surface with impurities and pre-existing 
broken inter- and intra-chain bonds, the simulated initial 
cellulose surface is perfectly regular and fully crystalline. 
Already, with only 5% of bonds hydrolyzed, the substrate 
becomes highly irregular (Shishir Chundawat, Personal 
Communication). In order to observe the effect of these 
irregularities, simulations were performed with an initial 
percentage of broken glycosidic bonds. Figure 9 shows 
how the initial slow hydrolysis phase diminishes as the 
initial cellulose substrate becomes more and more 
irregular. 

Even though we obtain reasonable qualitative agree- 
ment, it cannot be simply improved by just using a few 
experimentally determined kinetic parameters. Since 
our model is a detailed one, there are several other 
parameters that need to be optimized as well to get 
quantitative correspondence. For example, when we use 
the experimentally determined rates for adsorption 
(k on = 4.2 * 10 s" - based on association constant 
K a =1.4*10 6 M'V 1 [54] and k off [55]) and desorption 




Figure 9 Hydrolysis by endo-cellulases for imperfect crystals. 

Time course of hydrolysis by endo-cellulases when various glycosidic 
bond percentages are initially cleaved in the cellulose layer. The 
inset shows the same curves for a short time scale. 



(k off =0.03 M' 1 s" 1 [55]) for EG-I from cellulose crys- 
tals, we observe that the enzyme hydrolyzes the cellu- 
lose crystals completely in under 15 hours. Other 
physical reasons may also contribute. One of the rea- 
sons for the fast processing rate observed in the model 
is because we only process two-dimensional crystals of 
cellulose, while cellulose crystals are three-dimensional 
in nature. In three-dimensional crystals, there are mul- 
tiple layers of cellulose chains in the crystal, and not 
all of the glycosidic bonds are available as substrate for 
the enzyme to process from the beginning. Rather the 
inner layers of the crystal are available as substrate only 
after the layers above them are partially processed. 




20 30 40 

Time(h) 

Figure 10 Time course of hydrolysis by exo-R cellulases. The 

adsorption rate constant was set to 10 3 1/(sM)). The inset shows the 
variation in the conversion time needed to degrade 5%, 50% and 
80% of the substrate as a function of the adsorption rate constant 
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Another reason could be that the rate constant for 
decrystallization (or chain separation) of cellulose crys- 
tals by the enzymes is based on the theoretical time- 
scales for breaking of hydrogen bonds in an aqueous 
environment, which is quite rapid. 

Hydrolysis by exo-cellulases 

The hydrolysis profile produced by exo-cellulases is 
shown in Figure 10. The cellulose conversion is constant 
over time, a feature that has not been observed in 
experiments [16,50]. As expected, the time it takes to 
convert a given percent of the substrate decreases as k on 
increases (Figure 10, inset). 

The effect of varying the initial [E] 0 enzyme concentra- 
tion upon conversion times is plotted in Figure 11a. In 
contrast to Figure 8a, the substrate becomes saturated 
by exo-cellulases at lower enzyme concentrations than 
we observed in the case of endo-cellulases. This is be- 
cause the number of free chain ends always remains 
small compared to the number of enzyme particles in 
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Figure 1 1 (a) Effect of initial concentration of exo-cellulases. 

Change in convencion time as initial enzyme concentration is varied 
(b) Sugar production over time. System parameters are: N = 25000 
glucose units, k on = 10 5 1/(sM), k off = 100 1/s. 



solution. Figure lib shows that the only sugar produced 
by exo-cellulases is cellobiose. Experimental results, 
however, report the production of both glucose and cel- 
lotriose along with cellobiose, although cellobiose is the 
major product [16,50]. The relatively constant processive 
speed of exo-R cellulases along the glucan chain explains 
the constant hydrolysis rate observed in Figure 10 and 
the non-decreasing gap between the two curves in 
Figure 11a. 

The relatively constant processing time of the exo- 
cellulases is the result of using a simple coarse-grained 
description of processing time for exo-cellaloses since 
no rates have been measured for specific events asso- 
ciated with processivity. Values in these calculations are 
set such that k on , k off >> l/t move for exo-cellulases. Thus 
the binding of the enzyme to the substrate is at equilib- 
rium. These parameters can be easily adjusted to match 
up with any forthcoming experimental observations. In 
addition, the total concentration of the substrate (redu- 
cing ends for exo-R or non-reducing ends for exo-N) 
does not change with time until the whole chain is pro- 
cessed. This ensures that the concentration of enzyme- 
bound substrate remains nearly constant with time for 
the exo-cellulases resulting in a nearly constant rate of 
processing until the end. We expect these effects to be 
reduced in three-dimensional crystals of cellulose in 
which multiple layers of cellulose chains have to be pro- 
cessed by the exo-cellulases as the chains get hydrolysed 
in a staggered fashion in these crystals. It has also been 
reported [56] that cellulose-binding modules bind to in- 
soluble non-crystalline cellulose with a 10-20-fold 
greater affinity than to cello-oligosaccharides and/or sol- 
uble polysaccharides. Future expansion of this model 
will incorporate a non-constant adsorption rate of 
enzymes that would depend on the length of the cello- 
oligosaccharides; this will bring further complexity to 
the model. In addition, incorporation of stochasticity 
in the processing of the cellulose chain by the exo- 
cellulases and better estimates for rates of decrystalli- 
zation of the cellulose crystal could lead to better 
agreement with experimental hydrolysis rates. 

Hydrolysis by endo- and exo-cellulases 

In order to test whether our model reproduces the ex- 
perimentally observed endo-exo synergy, we used ex- 
perimental data reported by Eriksson and colleagues (see 
Figure 1A [50]) and modified the kinetic rate constants 
for both types of cellulases by fitting the model single 
enzyme hydrolysis curves to the experimental data. The 
initial enzyme concentration was set to [E] 0 = 1.5 fcM 
and the cellulose concentration to 10 g/L [50]. Numer- 
ical results show the endo-exo synergy (Figure 12a) and 
the time scale is the same as observed experimentally 
[50]. It should be kept in mind that the substrate in that 
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Figure 12 Endo-exo synergy, (a) Comparison of the time course of 
sugar production when cellulose is degraded only by endo-cellulases 
(red curve), only by exo-R cellulases (green curve), or when the two 
types of enzymes were acting together (orange curve). The sum of 
the red and green curves (blue curve) yields lower sugar production 
than in the case when the substrate is converted simultaneously by 
endo- and exo-cellulases. (b) Adsorbed enzyme percentages during 
hydrolysis when the cellulases are acting alone (red and green dots) 
or in a mixture (maroon and orange dots). 



experimental work is steam-pretreated spruce (lignocel- 
lulose), not pure cellulose. However, it is encouraging 
that we do observe similar behavior. Sugar production is 
higher when exo-R cellulases hydrolyze the cellulose sur- 
face in the presence of endo-cellulases when compared 
to the sum of their conversions achieved alone. 

The increase of free chain ends produced by endo- 
cellulases is the primary source of the endo-exo synergy 
observed in the model. Figure 12b illustrates how this ef- 
fect contributes to a large increase in the percentage of 
adsorbed exo-cellulases. This percentage is constant 
when exo-cellulases degrade the substrate alone, while 
in the presence of endo-cellulases it grows to higher 
values, contributing to a fast and efficient degradation of 
the substrate. 




ExoR 11:1 5:1 3:1 2:1 7:5 1:1 5:7 1:2 1:3 1:5 1:11 Endo 

ExoR/Endo 
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KxoR 5:1 2:1 1:1 1:2 1:5 Endo 

ExoR/Endo 

Figure 13 Effect of composition of endo-exo mixture, (a) k on 

(endo) = W(sM)' 1 , k off (endo) = 0.0} s'' and (b) k on (endo) = 100(sM)'\ k off 
(endo) = 0.1 s . In both cases 

k on (exo-R) = W 4 (sM)', k off (exo-R) = Ws' 1 . {N = 25000 glucose units, total 
cellulase concentration is 2 fjM). 



Results regarding synergism between pure Tricho- 
derma cellulases [14,57] showed that the endo-exo syn- 
ergy depends on the ratio of the concentrations of the 
individual enzymes. Here, we tested whether our model 
qualitatively reproduces this observation by comparing 
conversion times — the time to degrade 5%, 25%, 50% or 
80% of the substrate — for various exo-R/endo ratios. 
Using the hydrolysis rates listed in Table 4, we consider 
two cases: i) the overall hydrolysis of cellulose by endo- 
cellulases takes place at a slower rate than the overall 
hydrolysis of cellulose by exo-cellulases (Figure 13a); ii) 
the overall hydrolysis by exo-cellulases is set to be 
slower than that by endo-cellulases (Figure 13b). As the 
rate-limiting step in the model is the adsorption of cellu- 
lases onto the substrate, we attain this by varying the k on 
adsorption rate constant while fixing the equilibrium 
constant of each of the cellulases. In both cases the sub- 
strate conversion time has a minimum at a specific ratio 
of the concentrations of the individual enzymes. For the 
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first case (see Figure 13a) the optimal exo-R/endo ratio 
is 2:1, while for the second case (see Figure 13b) this 
ratio is 5:1. These minimum conversion times in both 
cases are much smaller than the conversion times 
obtained in single cellulase runs. These optimal ratios 
were obtained for a perfectly regular cellulose substrate; 
however, as pointed out earlier [14,57], the optimal ex- 
perimental ratio is strongly dependent on the character- 
istics of the substrate. For example, on filter paper [57] 
the optimal exo-R/endo ratio was found to be 74:26 
when the total enzyme concentration was 1 fiM and 
90:10 when the total enzyme concentration was 10 fj.M. 



Sensitivity analysis 

We have also carried out a sensitivity analysis to evaluate 
the relative importance of some of the physical quan- 
tities involved in the simulations. First, we studied the 
volume dependence of the overall hydrolysis. The reac- 
tion volume, V, determines the adsorption rate and 
therefore it determines the amount of adsorbed cellu- 
lases, affecting the rate of the overall hydrolysis. The 
simulation behaves as expected: a smaller volume results 
in higher adsorption rates, leading to faster hydrolysis 
(Figure 14). 

Next, the absolute size of cellulose surface was varied 
to determine the effect of the size of the simulation cell 
(Figure 15a). It is reassuring to see that the size of the 
substrate does not have any effect on the oligomer distri- 
bution. The amount of glucose, cellobiose and cellotriose 
increases as the substrate becomes larger, but their rela- 
tive concentrations are not affected by the substrate. A 
larger substrate needs more time to be degraded by the 
same amount of cellulases. The amount of adsorbed 
endo-cellulases is plotted as function of time in 
Figure 15b. The system size does not qualitatively 
change the number of adsorbed cellulases, it only affects 
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Figure 14 Volume dependence of the time course of hydrolysis 
by endo-cellulases [N E = 10). 



the overall hydrolysis time, and as such the curves get 
shifted towards larger time scales. 

Conclusions 

In an effort to complement both all-atom molecular dy- 
namics and coarse-grained simulation tools, we have 
developed a an agent-based the dynamical, functional 
model capturing the surface chemical reaction of cellulose 
hydrolysis by enzymes at the molecular level. This model 
accounts for heterogeneous enzymatic hydrolysis reactions 
occurring on the substrate surface (a reaction taking place 
in dimensions less than three), and incorporates key fac- 
tors controlling it that are different from those in an aque- 
ous environment. The catalysis process is broken down 
into distinct parts related to different kinetic events car- 
ried out by individual particles. These events are essen- 
tially chemical reactions taking place on the surface of 
cellulose (adsorption, breaking inter-chain hydrogen 
bonds, cleaving glycosidic bonds, desorption) and consti- 
tute the main elements of this model. Reactions are moni- 
tored by following and updating the state (based on a set 
of predefined rules) of each individual particle in the 
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Figure 15 The effect of varying the substrate size, (a) Time 
course of sugar production and (b) of adsorbed cellulases when 
cellulose is hydrolyzed by endo-cellulases (N £ = 10). 
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system. Simulation results showed good qualitative agree- 
ment with experimental data. 

The agreement with experimental data can be improved 
by obtaining better experimental estimates of the para- 
meters in Table 4 and by extending the current model to 
three dimensions. Initial experiments that can greatly 
benefit the model are those that probe the kinetics of dif- 
ferent steps for the individual domains of the cellulases 
(Carbohydrate Binding Domain and Catalytic Domain) 
separately. These experiments need to quantify the bind- 
ing affinities, k on and k 0 jf. Then similar measurements on 
the entire cellulases for binding of the same substrate will 
help to verify the role of individual domains and provide a 
measure of productive and non-productive binding. These 
measurements need to be carried out for pure cellulose 
substrates of different shapes morphology, degree of 
polymerization (DP) and partially digested states. 

Finally, our model only simulates the degradation of a 
single cellulose crystal layer, a feature that should be 
extended to capture the degradation of a whole cellulose 
crystal. The major effect from a three-dimensional 
model is expected to be that the substrate would not be 
completely accessible at the same time for the cellulases 
to digest. Also, such a three-dimensional model can cap- 
ture the possibility that floating sheets of detached sub- 
strate may slow cellulases from reaching a larger surface 
where more efficient digestion is possible. The current 
model does not account for surface diffusion, which is 
likely to be important based on results reported by Jervis 
and colleagues [53], who showed that diffusion does not 
limit enzyme activity. Fortunately, these deficiencies are 
not of a fundamental nature because our model is easily 
extendable and can incorporate them as well as add- 
itional properties of various cellulase systems on differ- 
ent types of cellulose surfaces. Importantly, this 
approach could be broadened to other classes of cellu- 
lases and even to cellulosomes as additional experimen- 
tal data becomes available. For this reason we believe 
that this model constitutes a significant contribution to 
the ability to simulate the complicated reactions 
involved in cellulose degradation. 

Endnotes 

a Each row corresponds to one of the seven parameters 
characterizing one monomer while the columns repre- 
sent numerical values the parameters can take. Each 
entry of the table denotes a distinct condition of a glu- 
cose unit. For a detailed explanation, please see text. 
b See reference [50]. c See reference [22]. d See reference 
[49]. 
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